Universality classes in directed sandpile models 
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We perform large scale numerical simulations of a directed version of the two-state stochastic 
sandpile model. Numerical results show that this stochastic model defines a new universality class 
with respect to the Abelian directed sandpile. The physical origin of the different critical behavior 
has to be ascribed to the presence of multiple topplings in the stochastic model. These results provide 
new insights onto the long debated question of universality in abelian and stochastic sandpiles. 
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The class of sandpile models, consisting of the original 
Bak, Tang and Wiesenfeld (BTW) 0] automata and its 
theme variations, is considered the prototypical exam- 
ple of a special class of driven non-equilibrium systems 
exhibiting a behavior dubbed self-organized criticality 
(SOC). Under an external drive, these systems sponta- 
neously evolve into a stationary state. In the limit of 
infinitesimal driving the stationary state shows a singu- 
lar response function associated to an avalanche-like dy- 
namics, indicative of a critical behavior. Sandpile models 
have thus attracted a great deal of interest, as plausible 
candidates to explain the avalanche behavior empirically 
observed in a large number of natural phenomena ||^. 

In recent years, the possibility of understanding the 
sandpile critical behavior in analogy with other non- 
equilibrium critical phenomena such as branching pro- 
cesses , interface depinning models |^,^ , and absorb- 
ing phase transitions (7[g] has been pointed out. It is then 
most important to identify precisely, for sandpiles, the 
universality classes and upper critical dimensions, which 
are basic and discriminating features of the critical be- 
havior. Despite significant numerical efforts, however, 
these issues remain largely unresolved. For instance, it is 
still an open problem wheter or not the original determin- 
istic BTW sandpile and the stochastic Manna two-state 
model belong to the same universality class. Theo- 
retical approaches ||lO|-p^ support the idea of a single 
universality class, while numerical simulations provide 
contradictory results |p^|-|l5|| . 

In order to have a deeper understanding of the uni- 
versality classes puzzle, we turn our attention to directed 
sandpile models [ p^ . In this case Dhar and Ramaswamy 
obtained an exact solution for the Abelian directed sand- 
pile (ADS) , that can be used as a benchmark to check 
the numerical simulation analysis. Directed sandpiles 
thus become an interesting test field to study how criti- 
cal behavior is affected by the introduction of stochastic 
elements. Despite the fact that results obtained for di- 
rected models cannot be exported "tout court" to the 
isotropic ones, the eventual apparence of different uni- 
versality classes provides interesting clues on the general 
problem of universality in sandpiles. This issue has been 
recently addressed in a particular case by Tadic and Dhar 
[O , but a general discussion of universality classes in di- 



rected sandpile automata is still lacking. 

In this letter we present large scale numerical simula- 
tions of Abelian and stochastic directed sandpile (SDS) 
models. First, we study an ADS model for which we re- 
cover numerically the results expected from the analyti- 
cal solution [|l6|. Then we introduce a stochastic model 
which is a directed version of the Manna two-state sand- 
pile In this case, the set of critical exponents defines 
a different universality class. For both models we provide 
a very accurate study of finite size effects and the conver- 
gence to the asymptotic behavior. For small and medium 
lattice sizes we find scaling anomalies that are similar to 
those encountered in isotropic models. We also study in 
detail the geometrical structure of avalanches. The pres- 
ence of multiple topplings appears to be the fundamental 
difference between Abelian and stochastic models. Nu- 
merical simulations in Euclidean dimension d > 2 show 
that both universality classes have an upper critical di- 
mension dc = 3, where strong logarithmic corrections to 
scaling are present. 

We consider the following definition for an ADS model 
(see Fig. |l|(a)): On each site of a c?-dimensional hyper- 
cubic lattice of size L, we assign an integer variable z^, 
called "energy" . Each time step, an energy grain is added 
to a randomly chosen site (z^ Zi + 1). When a site ac- 
quires an energy greater than or equal to the threshold 
Zc — 2d — 1, it topples. Topplings are directed along 
a fixed direction (defined usually as "downwards"): 
when a site on the hyperplane oijj topples, it sends deter- 
ministically one energy grain to each nearest and next- 
nearest neighbor site on the hyperplane X|| -I- 1, for a 
total of 2d — 1 grains. This definition differs from the 
ADS studied by Dhar and Ramaswamy in the ori- 
entation of the lattice. Both models, however, share the 
same universality class, being Abelian, deterministic, and 
directed. 

The stochastic generalization of the above model is de- 
picted in Fig. 0(b): The threshold is now Zc = 2, inde- 
pendent of the spatial dimension. When a site at the hy- 
perplane x\\ topples, it sends two grains of energy to two 
sites, randomly chosen among its 2d — 1 neighbors in the 
hyperplane -I- 1. The dynamical rule of this model can 
be defined exclusive if the two energy grains are always 
distributed on different sites. On the contrary, a nonex- 
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FIG. 1. Toppling rules in d = 2 for directed sandpiles. 
Filled circles represent active (toppling) sites; empty circles 
are stable sites. In the Abelian model (a) an active site sends 
one grain to each one of its three neighbors on the next down- 
wards row. In the stochastic models (b) one grain is sent to 
two randomly chosen downwards neighbors. 

elusive dynamics allows the transfer of two energy grains 
to the same site. We will consider separately the cases 
of the exclusive stochastic directed sandpile (ESDS) and 
the nonexclusive stochastic directed sandpile (NESDS). 
It is worth remarking that stochasticity does not alter the 
Abelian nature of the model [|l8|. All three models are 
locally conservative; no energy grains are lost during a 
toppling event. Boundary conditions are periodic in the 
transverse directions and open at the bottom hyperplane 
X|| = L, from which energy can leave the system. 

In the critical stationary state, we can define the prob- 
ability that the addition of a single grain is followed by an 
avalanche of toppling events. Avalanches are then char- 
acterized by the number of topplings s, and the duration 
t. According to the standard finite size scaling (FSS) hy- 
potesis, the probability distributions of these quantities 
are described by the scaling functions 



P{t) = t-^*Tit/t,) , 



(1) 
(2) 



where Sc and tc are the cutoff characteristic size and time, 
respectively. In the critical state the lattice size L is the 
only characteristic length present in the system. Ap- 
proaching the thermodynamic limit (L — > oo), the char- 
acteristic avalanche size and time diverge as Sc ~ and 
tc ^ L'^, respectively. The exponent D defines the frac- 
tal dimension of the avalanche cluster and z is the usual 
dynamic critical exponent. The directed nature of the 
model introduces a drastic simplification, since it imposes 
z = 1. A general result concerns the average avalanche 
size (s), that also scales linearly with L |l6| , p^ , p0[ : a new 
injected grain of energy has to travel, on average, a dis- 
tance of order L before reaching the boundary. In the 
stationary state, to each energy grain input must cor- 
respond, on average, an energy grain flowing out of the 
system. This implies that the average avalanche size cor- 
responds to the number of topplings needed for a grain 
to reach the boundary; i.e. (s) ^ L. The same result 
can be exactly obtained by inspecting the conservation 
symmetry of the model p^] . 



For the ADS, the exact analytical solution in d = 2 
yields the exponents = 4/3, Tt = 3/2 and D = 3/2 
]l6[ . The upper critical dimension is found to be dc = 3, 
and it is also possible to find exactly the logarithmic cor- 
rections to scaling ]l6| , p^ . The introduction of stochastic 
ingredients in the toppling dynamics of directed sand- 
piles has been studied only recently in a model that ran- 
domly stores energy on each toppling [p^ . This model is 
strictly related to directed percolation and defines a uni- 
versality class "per se" . In our case stochasticity affects 
only the partition of energy during topplings, and there 
is no analytical insight for the critical behavior of this 
model. In order to discriminate between ADS and SDS 
we perform simulations of both models for sizes ranging 
from L = 100 to L = 6400. Statistical distributions are 
obtained averaging over 10^ avalanches. Comparison of 
numerical results on the ADS allows us to check the re- 
liability and degree of convergence with respect to the 
lattice sizes used. 

It is well known from the many numerical papers on 
sandpiles that an accurate determination of the expo- 
nents Ts and Tt is a subtle issue. An overall determina- 
tion within a 10% of accuracy is a relatively easy task. 
However, a truly accurate measurement, allowing a pre- 
cise discrimination of universality classes, is strongly af- 
fected by the lower and upper cut-offs in the distribution. 
Extrapolations and local slope analysis are often very 
complicated and the relative error bars are not clearly 
defined. In this respect, it is far better to calculate ex- 
ponents by methods that contain the system-size depen- 
dence explicitly; namely data collapse and moment anal- 
ysis. Moment analysis was introduced by De Menech et 
al. in the context of the two dimensional BTW, and 
it has been used extensively on Abelian and stochastic 
models |l^j2^ . The g-moment of the avalanche size dis- 



tribution on a lattice of size i, 
the following size dependence 



= / s«P(s)ds, has 



where we have used the transformation y = s/L^ in 
the finite size scaling (FSS) form Eq. (0). More gener- 
ally, (s^)^ ~ L°'=('?\ where the exponents (Ts{q) can be 
obtained as the slope of the log- log plot of (s')^ ver- 
sus L. Using Eq. (|), we obtain (s«+i)i/(s«)L ~ or 
as{q+l) — crs{q) = D, so that the slope of (Ts{q) as a func- 
tion of q is the cutoff exponent; i.e. D — d(Js{q)/dq. This 
is not true for small q because the integral in Eq. (^) is 
dominated by its lower cutoff. In particular, corrections 
to scaling are important for q < Tg — l. An additional and 
strong check on the numerical data can be found in the 
fact that, as we have previously shown, the first moment 
of the size distribution must scale linearly with L. This 
last constrain also allows the evaluation of the the expo- 
nent Tg from the scaling relation (2 — Ts)D = (Ts(1) = 1, 
that should be satisfied for large enough sizes. 

Along the same lines we can obtain the moments of the 
avalanche time distribution. In this case (i'^)L ~ L"'*^'^\ 
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FIG. 2. Plot of as{q) for the d = 2 models ADS, ESDS, 
and NESDS. 



with dat{q)/dq = z. Analogous considerations for small 
q apply also for the time moment analysis. Here, an 
estimate of the asymptotic convergence of the numerical 
results is provided by the constraint z = 1, that holds for 
large enough sizes. Then, the Tt exponent can be found 
using the scaling relation (2 — Tt) = crt(l). 

Despite the fact that the moment method is usually 
rather accurate, it must be corroborated by a data col- 
lapse analysis. The FSS of Eqs. has to be verified 
and must be consistent with the numerical exponents ob- 
tained from the moment analysis. This can be done by 
rescaling s s/L^ and P{s) P{s)L^'^' and corre- 
spondingly t tjU and P{t) -> P{t)L^''K Data for 
different L must then collapse onto the same universal 
curve if the FSS hypothesis is satisfied. Complete consis- 
tency between the methods gives the best collapse with 
the exponents obtained by the moments analysis. In Ta- 
ble I wc report the exponents found for the ADS, ESDS 
and NESDS in d = 2. Figure ^ shows the moments (Ts(g). 
Figures || and ^ plot the FSS data collapse for sizes and 
times, respectively. 

The exponents obtained for the ADS are in perfect 
agreement with the expected analytical results. This 
fact supports the idea that the system sizes used in the 
present work allow to recover the correct asymptotic be- 
havior. It is worth remarking that, for small and medium 
lattice sizes, both moments and data collapse analysis 
present scaling features that can not be reconciled in the 
single scaling picture usually considered. These anoma- 



Model 


Ts 


D 


Tt 


z 


DR 
ADS 
ESDS 
NESDS 


4/3 
1.34 ±0.01 
1.43 ±0.01 
1.43 ±0.01 


3/2 
1.51 ±0.01 

1.74 ±0.01 

1.75 ±0.01 


3/2 
1.51 ± 0.02 
1.71 ±0.03 
1.74 ±0.04 


1 

1.00 ±0.01 
0.99 ±0.01 
0.99 ±0.01 


TABLE I. Critical 


exponents 


for directed 


sandpiles in 



CO 



O 



CO 



O 









a) 




1.34 








1.51 






6 


-4 

lo^ 


ho s/L'' 





- 
- 




1 ' 


1 ' 

b) - 




1.43 






D = 


1.74 


1 





logio s/L 



D 



d = 2. DR: Dhar and Ramaswamy's exact result; 
Abelian model; ESDS, NEDSD: stochastic models. 
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FIG. 3. Data collapse analysis of the integrated avalanche 
size distribution for the d = 2 models a) ADS and b) ESDS. 
System sizes are L = 400, 800, 1600, 3200, and 6400. The 
results for the ESDS and NESDS are identical, within error 
bars. 



lies are not persistent and disappear for reasonably large 
sizes {L ~ 10"^). This evidence for a slow decaying of 
finite size effects could shed light into several anomalies 
reported in isotropic sandpiles, for which, unfortunately, 
it is very difhcult to reach very large sizes [p^ , p3|j24|] . Re- 
sults for the ESDS and NESDS are identical within the 
error bars, indicating that these two models are in the 
same universality class. On the other hand, the obtained 
exponents show, beyond any doubt, that Abelian and 
stochastic directed sandpile models do not belong to the 
same universality class. 

The compelling numerical evidence for two distinct 
universality classes does not tell us what is the basic 
mechanism at the origin of the different critical behav- 
ior. In order to have a deeper insight into the dynamics 
of the various models, we have inspected the geomet- 
ric structure of the resulting avalanches. In Fig. |^ we 
depict in a color plot the local density of topplings in 
two avalanches of size 50 000 corresponding to the two- 
dimensional ADS and ESDS models. From the figure 
it becomes apparent that the stochastic dynamics intro- 
duces multiple toppling events, which are by definition 
absent in the Abelian case. This gives rise to very dif- 
ferent avalanche structures, eventually reflected in the 
asymptotic critical behavior. In particular, the fractal di- 
mension D is indicative of the scaling of toppling events 
with sizes. In the stochastic case we recover a higher 
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FIG. 4. Data collapse analysis of the integrated avalanche 
time distribution for the d = 2 models a) ADS and b) ESDS. 
System sizes are L = 400, 800, 1600, 3200, and 6400. The 
results for the ESDS and NESDS are identical, within error 
bars. 



fractal dimension than in the Abelian case. The multi- 
ple toppling mechanism has been proposed in the past 
as the origin of differences between isotropic Abelian and 
stochastic sandpiles as well. In that case, however, mul- 
tiple toppling is a common feature of both models, and 
for the largest sizes reached so far share, they share the 
same fractal dimension D |l5|| . 

Analysis of the models in three dimensions is strongly 
hindered by the presence of logarithmic corrections 
1^,^. Nonetheless, a naive application of the moment 
analysis yields values compatible with the mean-field re- 
sults Ts = 3/2^ = 2, and D = 2 More interest- 
ingly, in Ref. jl^l the authors were able to deduce the 
exact form of the logarithmic corrections in d = 3 for 
the avalanche time distribution, namely P{t) ~ t~^\Tit. 
In Figure ^ we have checked that the same logarithmic 
corrections apply to both the Abelian and ESDS sand- 
piles. This remarkable fact lends support to the critical 
dimension of the stochastic model being dc — 

In summary, we have reported large scale numerical 
simulations of a stochastic directed sandpile model. This 
model defines unambigously a different universality class 
with respect to the Abelian directed sandpile model. The 
origin of this difference is traced back to the avalanche 
cluster geometric structure, providing new clues to un- 
derstand the effect of stochastic elements in the dynamics 
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FIG. 5. Color plots of the local density of topplings in two 
avalanches of size 50000, in d = 2, for a) ADS and b) ESDS. 
Yellow represents a single toppling per site; red stands for the 
maximum number of topplings. 



of avalanche processes. 
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